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1.  Introduction 


This  final  report  covers  1  year  of  ARO-sponsored  research  into  the 
fundamental  behavior  of  rotor  dynamics.  The  original  proposal  was  for  3 
years,  but  the  first  two  years  were  done  under  a  separate  grant.  Thus,  this 
report  covers  only  a  portion  of  the  originally-proposed  work.  In  this  final 
report,  we  will  summarize  the  work  done,  including  publications  and 
scientific  personnel;  and  we  will  provide  pertinent  technical  descriptions  of 
each  major  area  we  have  pursued. 

The  philosophy  of  our  research  has  been  to  increase  the  fundamental 
understanding  of  those  dynamic  and  aerodynamic  phenomena  associated  with  the 
helicopter.  Our  approach  has  been  to  follow  three  intertwining  lines  of 
research.  The  first  of  these  is  the  line  of  mathematical  modeling.  Here,  we 
wish  both  to  synthesize  and  refine  mathematical  models  for  various,  isolated 
rotorcraft  phenomena,  and  to  learn  to  couple  them  together  in  a  systematic 
way.  This  is  the  building-block  approach  we  have  followed.  The  second  line 
of  inquiry  has  been  in  the  development  of  solution  methodologies  for  these 
equations.  Here,  certain  solution  strategies  work  better  for  certain  models; 
and  some  modeling  techniques  require  new  solution  strategies.  We  look 
specifically  at  methods  that  magnify  our  insight,  are  computationally 
efficient,  and  that  can  be  extended  to  large-scale  systems. 

This  leads,  then,  to  the  third  thread  of  research:  basic  physical 
insight.  Of  course,  because  we  deal  with  isolated  components  or  with 
simplified  couplings,  we  do  not  intend  to  be  able  to  make  predictions  on 
helicopter  stability  and  response  that  would  be  applicable  to  detailed  design 
studies.  On  the  other  hand,  we  do  expect  our  methods  to  be  predictive  of  the 

behavior  of  simplified  research  models,  such  as  those  used  by  the  Army 

Research  and  Technology  Laboratories.  Furthermore,  we  believe  our  results 
give  qualitative  insight  into  the  physical  phenomena  present  in  production 
rotors.  Thus,  we  try  to  involve  all  three  elements  in  our  research  effort. 

2.  Statement  of  Problem 

The  objectives  and  scope  of  this  work  are  as  follows: 

1)  To  discover  the  basic  relationships  between  blade  structural 
parameters  and  the  flap-lag-torsion  airloads  that  result. 

2)  To  determine  the  extent  to  which  rotor-body  coupling  affects 
inplane  loads  and  overall  helicopter  vibrations. 

3)  To  develop  our  basic  trim  procedures  to  the  point  at  which  they  can 

be  applied  to  large,  state-of-the-art  rotor  response  program. 

4)  To  determine  the  effect  of  dynamic  stall  on  the  rotor  airloads  and 
on  the  basic  trimming  methods. 

5)  To  investigate  other  methods  of  obtaining  time  histories  of  rotor 
response,  including  Hamilton's  Law  of  Varying  Action. 
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Before  proceeding  to  the  details  of  each  objective,  it  is  informative  to 
outline  the  scope  in  each  task.  With  the  exception  of  item  3,  the  above 
objectives  are  not  aimed  at  the  quantitative  prediction  of  helicopter 
response.  They  are  aimed  at  obtaining  fundamental  insight  into  how  rotor 
vibrations  develop  and  into  how  they  can  be  efficiently  calculated.  Thus,  in 
item  1  we  consider  a  simple  elastic-blade  model  with  elastic  flap,  lag,  and 
torsion.  Although  other,  more  sophisticated  flap-lag-torsion  models 
certainly  exist,  they  have  not  been  obtained  under  the  same  assumptions  nor 
with  the  same  purpose  in  mind  as  ours.  Thus,  we  have  proceeded  slowly  and 
carefully  to  make  sure  we  understand  the  physical  processes  at  each  step. 

In  item  2,  we  are  looking  at  a  fuselage  with  5  rigid  body  modes  and  4 
elastic  modes  (as  in  our  prior  work)  but  with  a  more  detailed  rotor  model. 
Naturally,  a  true  fuselage  will  have  many  more  elastic  modes;  but  we  look  at 
a  generic  frequency  sweep  that  could  be  representative  of  several  potential 
modes.  Since  we  have  already  found  that  flapping  motions  drive  inplane 
motions  (while  inplane  effects  flapping  much  less)  we  make  several 
simplifying  assumptions  to  increase  the  productivity  (and  physical 
interpretations)  of  the  work. 

Item  3  is  the  only  area  in  which  we  approach  the  area  of  applications. 
These  trim  procedures  are  now  fairly  well  understood  in  terms  of  theory,  and 
the  advancements  now  come  through  more  sophisticated  applications. 

Therefore,  we  have  reformulated  the  trim  procedures. 

Item  4  is  a  new  area  of  research  that  developed  out  of  our  dynamic-stall 
work.  It  is  not  in  our  scope  to  develop  any  dramatically  new  dynamic  stall 
procedures.  We  merely  take  existing  methodologies,  investigate  how  they 
should  be  modified  to  be  useful  for  simplified  vibration  analyses,  and  study 
the  resultant  effects  on  the  types  of  calculations  we  are  making. 

Item  5  is  also  a  new  area  of  research  which  developed  out  of  our  prior 
trim  investigations.  For  nearly  linear  systems,  the  trim  method  of  periodic 
shooting  is  equivalent  to  finding  and  inverting  the  Floquet  transition 
matrix.  (An  earlier  solution  method  in  our  research  also  relied  on  Floquet 
theory  for  vibration  analysis.)  Thus,  it  is  natural  to  look  for  more 
efficient  means  of  finding  the  transition  matrix.  One  possibility  is  the  use 
of  Hamilton's  Law  of  Varying  Action  with  comparison  functions  in  time.  In 
this  research  we  study  Hamilton's  Law  in  detail  with  respect  to  convergence 
and  efficiency. 
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5.  Summary  of  Results 


In  this  section,  we  summarize  the  results  of  our  research  for  this 
one-year  effort.  However,  (the  first  two  years  were  performed  at  Washington 
University  in  St.  Louis),  we  will  also  refer  to  how  this  third  year  completes 
the  work  begun  in  the  first  two  years.  Reference  numbers  refer  to 
publication  in  Section  4. 

5.1  Finite  Elements  in  Time 


Recently,  much  attention  has  been  give  to  numerical  application  of 
Hamilton's  Law  of  varying  action.  Hamilton's  Law  is  a  variational  statement 
about  "action"  which  provides,  for  the  time  domain,  what  variation  of  work 
provides  in  the  space  domain.  Thus,  these  applications  of  Hamilton's  Law 
result  in  finite  elements  over  the  time  domain;  and  these  can  be  either 
p-version,  h-version,  or  a  combination  of  the  two  (depending  on  the  choice  of 
test  functions).  However,  numerical  applications  of  Hamilton's  Law  have 
sometime  resulted  in  solutions  that  do  not  converge  as  the  number  of  elements 
(i.e.,  polynomials)  is  increased.  In  this  research,  a  convergence  proof  was 
found,  based  on  the  bilinear  formulation,  which  demonstrates  that  some 
formulations  are  not  truly  bilinear  and  may  not  converge.  The  proof  also 
leads  to  an  alternate,  bilinear  formulation  of  Hamilton's  Law  for  which 
convergence  is  assured.  The  bilinear  formulation  also  leads  to  an 
alternative  statement  about  dynamics.  In  particular,  the  "virtual  action" 
plus  the  variation  of  action  over  a  space-domain  must  always  sum  to  zero. 

Numerical  application  of  the  correct  bilinear  formulation  leads  to 
Lagrange  multiplier  with  the  physical  connotation  of  an  end  momentum  (which 
is  the  analogy  of  end  force  in  spatial  problems).  Thus,  initial  velocity  is 
treated  as  a  "natural"  rather  than  as  "geometric"  boundary  condition;  and  the 
Lagrange  multiplier  converges  to  the  unknown  momentum  (i.e.,  velocity)  at  the 
end  of  the  time  period.  Thus,  the  bilinear  formulation  is  a  "mixed  method". 
Accuracies  of  solutions  with  the  Lagrange  multiplier  are  an  order  of 
magnitude  better  than  those  which  use  the  derivative  of  shape  functions  for 
velocity. 

In  the  limit  as  one  takes  many  elements  with  only  a  few  polynomials 
each,  this  formulation  reduces  to  a  classical  time-marching  method,  (an 
h-version  finite  element)  similar  to  Euler,  Runge-Kutta,  or  predictor 
correctors.  In  the  limit  as  many  polynomials  are  used  per  element,  but  with 
only  a  few  elements,  the  method  becomes  similar  to  a  Ritz-Galerkin  procedure 
in  time  (  a  p-version  finite  element).  Results  show  that,  for  any  given 
problem  (as  characterized  by  the  computational  cost  of  a  function 
evaluation),  there  is  an  optimum  choice  of  polynomial  number  in  order  to  meet 
any  error  criterion  with  minimum  computational  effort.  Similarly,  depending 
on  the  problem,  a  particular  choice  of  polynomial  number  may  or  may  not  be 
more  efficient  than  conventional  time-marching  methods.  In  general,  finite 
elements  in  time  become  more  efficient  than  marching  as  the  desired  accuracy 
becomes  exacting  and  as  function  evaluations  become  computationally 
expensive. 


The  details  of  this  work  can  be  found  in  References  5  and  7.  Figures 
1-5,  taken  from  Reference  5,  summarize  the  major  conclusions  of  the  work. 
Figure  1  shows  the  blade  flapping  angle,  p,  at  the  end  of  one  rotor 
revolution  as  calculated  both  by  Hamilton's  Finite  Elements  and  by  our  new, 
bilinear  formulation.  One  can  see  that,  for  Hamilton's  Law,  p  fails  to 
converge  uniformly  as  the  number  of  basis  functions  is  increased.  Figure  2 
further  shows  that  this  divergence  is  not  restricted  to  a  certain  advance 
ratio.  Results  with  Hamilton's  Law  can  be  in  error  over  100%  even  with  12 
basis  functions  at  p  =  0.5.  The  new  formulation,  on  the  other  hand, 
converges  in  all  cases. 

Figure  3  provides  hp-optimization  information.  The  figure  gives  log.„ 
of  the  error  as  a  function  of  the  number  of  floating-point  multiplication^ 
required  in  the  computation,  M.  The  straight-lines  are  for  various  number  of 
polynomials  per  element,  n,  with  the  step  size  being  a  running  parameter 
along  the  lines.  One  can  see  that,  for  any  given  error  tolerance,  there  is  a 
minimum  M  given  by  the  interior  of  optimum  curve  formed  by  the  locus  of 
straight  lines.  The  x's  are  results  from  Hammings  predictor-corrector.  For 
this  case  (C=16  implies  3  to  4  sine  or  cosine  functions  in  each  coefficient), 
the  hp  finite  elements  are  always  superior  to  conventional  time  marching. 

Another  interesting  result  of  our  research  can  be  seen  in  plots,  like 
Figure  4,  which  show  error  on  the  interior  of  a  large,  p-version  element. 

The  error  norm  is  zero  at  t=0  (due  to  known  initial  values).  It  goes  through 
some  erratic  oscillations  (due  to  truncation  errors  near  x=0),  but  settles 
down  over  the  rest  of  the  element.  The  minimum  error,  however,  is  found  at 
t=T,  the  end  of  the  period.  The  values  at  the  end  of  the  period  are  exactly 
what  is  needed  for  Floquet  theory.  Thus,  finite  elements  in  time  (when 
correctly  formulated)  are  ideal  for  Floquet  applications. 

Finally,  we  consider  numerical  stability  when  finite  elements  are  used 
to  march  indefinitely,  Figure  5.  The  values  p  and  o  are  the  system  frequency 
and  damping  multiplied  by  the  element  length,  At.  The  exterior  of  the  large 
semi-circles  (as  well  as  the  interior  of  the  small  semi-circles  on  the  p 
axis)  are  conditions  of  numerical  instability,  "n"  is  the  number  of 
polynomials  per  element.  To  put  this  in  perspective,  typical  radii  of 
convergence  for  other  methods  are  near  2  =  1.4,  smaller  than  even  the  n=2 
result.  At  n  =  6,  a  step  size  equal  to  two  periods  (r=3n)  is  required  before 
instability  occurs.  From  Figure  3,  however,  we  see  that  such  a  large  step 
size  would  result  in  very  large  errors  with  or  without  the  instability.  No 
optimal  point  on  Figure  3  is  unstable. 

5.2  Dynamic  Stall 


A  major  portion  of  our  research  has  dealt  with  the  introduction  of  a 
modified  ONERA  model  (for  dynamic  stall)  into  rotor  elastic-blade  analysis. 
In  this  research,  reported  in  References  4  and  6,  flap  and  inplane  bending 
are  described  by  two  nonlinear,  partial-differential  equations  which  are 
coupled  together  as  derived  in  Reference  1.  Each  equation  consists  of  lower 
and  higher  order  terms.  Approximate  solution  methods  are  applied  to  these 
aeroelastic  equations.  The  aerodynamic  circulations  are  expressed  by  our 
unified  theory.  Only  lower  order  terms  of  the  equations  in  Ref.  1  are 
considered  in  the  analyses.  This  allows  the  investigation  of  the  general 
trends  of  the  result  without  going  into  cumbersome  computation  with  little 
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improvement.  In  the  following,  we  write  the  system  of  equations  to  be  solved 
including  lower-order  terms  only 

Flap  equation,  order  e 


—  +  A6iv 


—  +  +  +  +  ,  — 


~(r ,  +  f2)(.t  +  n  sin  ip) 


(1) 


Lag  equation,  order  e2 


— 

_  ( - -■■  ■■■  P+)  +  A5V  +  (A5  —  A6)(#  +  <t>)w+  +  +  +  —  2  f  iv  +  w^iix  +  v  —  vT  = 

-  •'O 


-7— 1  +  r  1  )(w+  n  cos  ip  +  V '  +  \v)  —  -d—CD(x  +  fi  sin  t p)~ 
oa  6  a 


(2) 


f,  equation,  order  1 


k  r(  +  -if  |  =  kaUy  +  5b e+ 


(3) 


f:  equation,  order  1 


k'V  1  +  2dwkT2  +  iv-(l  +  c/'ifi  =  —  W(1  •+•  t/‘)Ct/TAC.  +  ek(UxAC.  + 


cAC. 

£0 


c/,)] 


(4) 


It  is  interesting  here  to  note  that  the  terms  of  the  lag  equation  are  of 
higher  order  compared  to  the  terms  of  the  flap  equation.  Thus,  to  first 
order,  flap  drives  lag  in  an  important  way,  but  lag  has  a  lesser  effect  on 
flap. 


An  effective  way  of  obtaining  accurate  pitch  inputs  is  based  on  the  use 
of  an  automatic  feed-back  system  that  can  trim  the  helicopter.  The 
auto-pilot  equations  are  developed  and  expressed  as  follows: 
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This  automatic-feed-back  system  provides  control  to  the  helicopter  for 
numerical  purposes.  It  adjusts  the  pitch  of  the  blade  to  maintain  thrust, 
roll  moment,  and  pitch  moment.  The  parameters  Ag  and  are  controller 

gains.  The  parameters  tq  and  t1  are  time  constants.  The  grouping  8  (P 


is  a  coupling  parameter  giving  the  pitch-roll  coupling  of  0s  and  ©c>  The 
2 

grouping  is  an  estimate  of  rotor  thrust  in  the  absence  of  an  explicit 

Cy  equation. 


In  the  computational  algorithm,  input  parameters  define  the  basic 
aircraft  configuration  and  flight  conditions.  The  blade  type  is  defined  by 
its  root  stiffness,  solidity.  Lock  number,  damping,  and  airfoil  type.  Flight 
conditions,  on  the  other  hand,  are  characterized  by  the  advance  ratio,  and 
thrust  coefficient.  In  Table  1,  we  present  the  common  parameters  used  in  the 
cases  discussed  in  this  report.  These  parameters  are  selected  from  current 
helicopter  data  and  used  for  illustration  purposes.  It  is  noted  that  we  use 
only  a  constant  value  of  CQ.  In  future  work,  we  intend  to  include  a  variable 
Cq  based  on  stall  assumptions.  Thus,  the  major  study  here  is  the  effect  of 

stall  on  lift  but  not  on  drag. 


Vibration  analysis  includes  flap  and  lag  responses  and  their  sensitivity 
to  advance  ratio  and  thrust  coefficient.  The  torsion  effect  is  neglected  in 
this  study.  The  variation  of  the  automatic  control  settings  required  for 
trim  is  determined  for  a  blade  revolution  in  a  steady  state.  In  addition, 
the  change  of  these  settings  as  a  response  to  blade  stall  is  identified. 

Last,  the  vibration  in  the  flap  direction  obtained  from  the  stall  model  used 
in  this  research  work  is  compared  to  results  obtained  when  a  linear 
aerodynamic  model  is  used. 
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Physical  Description 


time  constants 


In  order  to  determine  the  response  of  a  helicopter  blade,  one  must 
determine  the  pilot  controls  (collective  pitch,  cyclic  pitch,  and  shaft  tilt) 
that  are  applied  to  the  rotor.  There  are  two  primary  methods  of  doing  this. 
One  method  is  based  on  formulas  derived  by  a  harmonic  balance  of  the 
equations  without  stall.  In  the  other  method,  an  automatic  pilot  is 
implemented  to  update  the  controls,  equations  (5-7).  In  either  case,  the 
purpose  of  the  controller  is  to  eliminate  first  harmonic  moments  (^n  the 
rotating  frame)  which  come  through  (to  the  fixed  frame)  as  steady  pitching 
and  rolling  moments.  In  other  words,  "trim"  implies  (among  other  things)  the 
balancing  of  these  moments. 

Figures  6-9  provide  plots  of  the  flapping  response  at  x  =  0.7.  Because 
this  is  a  rigid  blade  with  root  spring,  flap  displacement,  w,  is  a  direct 
measure  of  root  moment.  The  four  figures  show  a  variety  of  thrust 
coefficients  and  advance  ratios,  some  with  stall  (high  y,  high  CT),  and 
others  with  little  or  no  stall  (low  y,  low  CT).  We  see  immediately  that  the 
approximate  formulas  for  trim  give  a  large  amount  of  1/rev  (i.e.,  sin  iji)  in 
the  response.  Thus,  they  are  not  accurate  due  to  the  aspects  of  the  model 
neglected  in  such  formulas  (inplane  motion,  unsteady  aerodynamics,  nonuniform 
inflow  distribution,  and  stall).  The  results  with  the  automatic  pilot, 
however,  show  two  peaks  per  cycle,  indicative  of  little  1/rev,  mostly  2/rev, 
and  some  3/rev,  which  is  indicative  of  trim. 

It  is  interesting  that  most  investigations  in  this  area  have  had  trouble 
obtaining  a  correct  trim.  For  example.  Figure  10  shows  results  at  y  =  0.4 
from  Friedmann.  Two  different  iteration  schemes  are  used  for  the  two  plots. 
Notice,  however,  that  the  results  (although  labelled  as  "propulsive  trim") 
have  a  large  1/rev  component  very  similar  to  that  seen  in  our 
approximate-formula  results  of  Figure  6.  This  is  to  be  expected  because  only 
an  approximate  formula  is  used.  Figure  11  presents  other  results  from 
Chopra.  Here,  we  see  a  very  large  fore-to-aft  1/rev  for  a  rotor  supposedly 
in  "propulsive  trim"  at  C  =  0.1,  y  =  0.2.  Furthermore,  the  curves  show  that 
the  response  is  not  even  jSeriodic,  as  the  slopes  do  not  match  at  i|>  =  0°  and 
360°.  In  Figure  12  taken  from  a  later  Chopra  reference,  the  authors  attempt 
to  correct  the  lack  of  trim  found  in  Figure  11.  The  dashed  curve  is  the  old 
result  (although  in  this  paper  it  is  corrected  to  be  periodic),  and  the  solid 
curve  is  the  new  result.  In  this  new  result,  the  trim  solution  is  modified 
to  include  elastic  twist.  Notice,  however,  that  although  the  1/rev  has  been 
reduced  by  about  50%,  it  is  still  very  much  present.  Furthermore,  the  solid 
curve  is  not  all  periodic,  with  an  error  of  over  100%  in  the  slope  between  <|> 

=  0°  and  i>  =  360°.  This  points  out  the  difficulties  in  finding  a  good  trim 
solution  with  elastic-blade  equations.  The  autopilot  aids  greatly  in  this 
regard. 

In  summary,  unified-aerodynamic  model  has  been  introduced  in  the 
elastic-blade  equations.  This  model  is  an  extension  of  the  0NERA  lift  model 
and  includes  plunge,  unsteady  free  stream,  and  large  angles  of  attack.  An 
ordering  scheme  has  been  used  to1  segregate  the  important  terms  from  others. 
The  elastic-blade  equations  are  presented  as  lower-order  terms  and 
higher-order  terms  in  the  flap,  lag,  and  circulation  equations.  A  solution 
method  based  on  a  modified  Galerkin's  method  is  used  to  separate  the  time  and 
space  variables  in  the  differential  equations.  A  numerical  solution  is 
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Figure  12  Flap-Lag  Response 


obtained  by  solving  the  equations  for  time  history.  Two  methods  have  been 
used  for  trimming  purposes.  Dynamic  response  has  been  conducted  for  a 
variety  of  thrust  coefficients  and  advance  ratios.  Results  lead  to  the 
following  conclusions: 

1)  The  approximate  method,  in  which  pitch  settings  are  approximated, 
gives  large  once-per-rev  oscillations  in  the  flap  response. 

2)  Automatic  controllers  eliminate  the  once-per-rev  and  are  used 
successfully  beyond  stall  limit. 

5.3  Tip-Loss  Aerodynamics 

In  this  area,  we  have  made  two  very  important  contributions  to  the 
understanding  of  tip  loss.  First,  we  have  developed  methods  to  optimize  the 
lifting-surface  mesh.  This  optimization  not  only  improves  accuracy  but  also 
guarantees  convergence,  something  that  cannot  be  said  of  conventional  mesh 
choices.  Second,  we  have  developed  closed-form  estimates  of  the  far-wake 
contributions  to  induced  velocity  even  with  wake  contraction.  This  improves 
convergence  on  velocity  computations  by  a  factor  of  10.  These  two  methods 
together  from  the  nucleus  of  our  tip-loss  computer  code  which  is  now  being 
utilized.  Because  neither  of  these  developments  has  yet  been  published 
(although  papers  are  now  in  review),  we  attach  to  this  report  an  extended 
abstract  which  outlines  the  developments  in  mesh  choice;  and,  under  separate 
cover,  we  have  sent  to  ARO  a  manuscript  on  the  far-wake  methodology,  which 
has  been  submitted  to  the  Journal  of  Aircraft. 

5.4  Rotor-Body  Coupling 

Our  research  in  rotor-body  coupling  has  been  severely  hampered  by  a 
tragic  accident  involving  Huang  Ming-Sheng,  the  graduate  research  assistant 
involved  in  this  work.  Reference  8,  attached  to  this  report  summarizes  the 
basic  theory  behind  the  work,  which  was  developed  by  Mr.  Huang  prior  to  his 
accident.  Since  recovering  from  these  injuries,  he  has  undertaken  the  task 
of  coding  this  theory.  That  work  is  now  completed.  Due  to  the  loss  of  time, 
we  are  not  now  able  to  deliver  extensive  numerical  calculations  as  we  had 
hoped.  We  will,  however,  continue  this  research  (hopefully  under  future  ARO 
funding) . 

5.5  Dynamic  Inflow 

We  have  made  considerable  progress  in  this  grant  toward  achieving  more 
widespread  use  of  the  theory  of  dynamic  inflow.  Reference  2  provides  the 
ultimate  correlation  between  our  ARO-sponsored  inflow  model  and  wind-tunnel 
test  data  in  the  Ames  7x10  wind  tunnel.  Reference  3  gives  a  review  of  the 
past  30  years  of  dynamic  inflow  modeling  including  the  present  state  of  the 
art  and  current  research  interests.  Current  research  interests  which  we  have 
pursued  in  this  grant  are:  1)  role  of  tip-path  dynamics  on  actuator-disk 
assumptions,  2)  effect  of  hub  motions  on  dynamic  inflow  and  3)  relationships 
between  higher-harmonic  dynamic  inflow  and  unsteady  aerodynamics.  In  the 
first  item,  we  now  believe  that  tip-path  dynamics  should  not  alter  the 
momentum  statement  of  dynamic-inflow  theory.  In  the  second  area,  hub  motions 
can  effect  the  formulation  of  1 inearized  dynamic  inflow.  We  have  worked  out 


the  details  in  Ref.  3.  In  item  3,  we  have  uncovered  some  very  strong 
relationships  between  Loewy  theory  and  dynamic-inflow  theory.  Those  may 
allow  us  to  develop  a  3-dimensional,  unified  theory  of  unsteady  aerodynamics 
that  is  applicable  to  rotary-wing  problems.  An  extended  abstract  concerning 
this  new  theory  has  been  submitted  to  the  AHS  Forum  and  has  been  sent  to  ARO 
under  separate  cover. 

5.6  Multiblade  Transform 

One  of  the  smaller  tasks  in  our  research  effort  has  been  the  study  of  a 
modified  multiblade  transform  designed  to  put  the  differential  mode  into  the 
nonrotating  system.  This  work  is  essentially  complete,  although  the  student 
involved  has  not  yet  finished  writing  his  M.S.  thesis.  The  major  conclusions 
are  as  follows: 

1)  The  new  transform  for  2-bladed  rotors  does  provide  improvement  over 
the  conventional  2-bladed  transform.  In  particular,  it  captures  the 
essentials  of  the  1/rev  instability  at  high  y. 

2)  For  rotors  with  an  even  number  of  blades  greater  than  2(e.g., 
4,6,8,etc.),  there  are  two  possible  alternate  transforms.  One  provides  some 
improvement  in  the  transient  analysis  but  degrades  accuracy  of  the  forced 
response.  The  other  improves  forced  response  but  degrades  transients. 
Therefore,  it  is  doubtful  that  either  transform  will  be  of  general  use  in 
dynamics  analyses. 

3)  As  a  spin-off  of  the  study,  we  have  developed  matrix-manipulation 
algorithms  that  can  perform  the  multi-blade  transformation  (either 
conventionally  or  in  one  of  our  new  ways).  These  algorithms  offer  a  simple 
way  to  make  this  important  change  of  variable,  and  they  are  superior  to 
computer  algebra,  numerical  analysis,  or  fourth-order  tensor  approaches  used 
by  other  investigators. 


Summary  of  Results 


In  summary,  we  can  say  that  our  research  has  been  very  successful.  It 
has  resulted  in  improved  analysis  tools  and  in  improved  understanding  of 
rotor  vibrations.  Many  of  these  tools  have  been  (and  are  being)  integrated 
into  production  rotor  codes. 
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ABSTRACT 


I 


Researchers  often  use  lifting-line  and  lifting-surface  (or  panel)  theories  to  obtain 
lift  and  drag  of  fixed  and  rotating  wings.  The  choice  of  panels,  line  segments,  and  control 
points  within  panels  (i.e.,  collation  points)  is  usually  performed  on  an  ad  hoc  basis  based  on 
engineering  judgement  The  results  of  this  research  show  that  often-used  methods  such  as:  1) 
placing  control  points  at  mid-span  of  the  panel,  2)  placing  large  panels  near  the  blade  root 
and  small  panels  near  the  wing  tip.  and  3)  using  equally-spaced  panels  near  the  blade  tip,  all 
result  in  order-unity  errors  in  the  calculation;  and  these  errors  do  not  decrease  with  refined 
mesh  even  as  the  number  of  panels  and  the  number  of  significant  digits  is  increased  without 
bound.  They  are  non-vanishing  residuals. 

The  work  reported  here  describes  both  how  to  eliminate  these  residual  errors 
through  proper  choice  of  panel  size  and  control  points,  and  how  to  optimize  the  mesh  size 
to  give  the  maximum  rate  of  convergence  as  the  number  of  panels  is  increased.  Comparisons 
with  experimental  data,  with  other  lifting-line  and  lifting-surface  results,  and  with  a 
classical  Fourier  solution,  demonstrate  the  superiority  of  the  new  procedure  over 
conventional  mesh-choice  methods. 


INTRODUCTION 


It  is  commonly  agreeed  that  the  key  to  accurate  calculation  of  the  rotor  aerodynamic 
behavior  is  the  correct  modeling  of  the  rotor  wake.  McCroskey1  concluded  that  lifting-line 
calculations  (for  rotary-wings)  are  in  error,  regardless  of  the  wake  model.  The 
prescribed-wake  method  of  calculation  also  gives  errors  in  spanwise  loading  distribution, 
regardless  of  which  representation  is  used  for  the  surface  of  the  blade.  Even  with  the 
complicated  free-wake  models,  lifting-surface  codes  still  fail  to  predict  adequately  some 
cases  with  highly  nonlinear  twist  distribution.  It  has  been  proposed  that  a  major  part  of 
these  observed  discrepancies  is  due  to  improper  selection  of  collocation  points,  (i.e.,  the 
possibility  of  running  into  mathematical  singularities  that  have  no  physical  counterpart2). 

Most  present-day  rotor  analyses  employ  a  Kernel-function  (in  the  Mangier  sense), 
which  contains  a  higher-order  singularity  and  is  diffcult  to  handle.  Some  researchers, 
however,  employ  the  Vortex  Lattice  Method  (VLM)  which  amounts  to  a  Cauchy-type 
finite-element  solution  to  both  lifting-surface  and  lifting-line  problems.  Many  applications 
of  the  VLM  have  been  made  to  problems  of  the  aerodynamic  analysis  and  design  of  wings 
with  considerable  sucess.3'*  However,  the  prediction  of  the  detailed  aerodynamic 
performance  of  a  rotary  wing  is  more  difficult  than  that  of  a  fixed  wing.  In  the  latter,  the 
wake  trails  back  from  the  wing  in  a  relatively  straight  path  to  downstream  infinity.  The 
effect  of  the  trailing  vortex  on  the  calculation  can  be  minimized  by  use  of  certain 
mathematical  techniques  or  by  use  of  a  free-wake  model,  which  automatically  allows  wake 
roll-up.  In  the  case  of  a  rotary  wing,  the  blades  pass  directly  over  their  own  wakes  as  well 
as  those  of  other  blades  as  they  rotate.  Furthermore,  a  given  element  has  a  longer  residence 
time  in  the  immediate  vicinity  of  the  rotor  plane  compared  with  that  of  a  fixed-wing  vortex 
element.  Many  computational  efforts  has  been  performed  to  reduce  the  error  caused  by 
wake  vortex  elements  of  rotary  wings.  These  include  the  numerical  integration  technique' 
(for  the  prescribed-wake  models),  curved-vortex  elements  8  (for  free-wake  models),  and 
division  of  the  wake  into  three  separate  regions  with  each  computed  separately.9 
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The  intention  of  this  research  is  to  discover  if  a  simple,  systematic,  optimized  VLM 
approach  can  be  developed  for  application  to  both  fixed  and  rotary  wings  (especially  for 
helicopter  blades),  despite  criticisms  of  the  VLM.  which  continue.  These  criticisms  contend 
either  that  the  lattics  are  laid  out  in  a  preconceived  manner  to  give  some  desired  answer  or 
that  too  many  lattices  are  required  for  adequate  convergence  of  the  computed  leading.  The 
present  work  is  to  derive  systematically  an  optimized  vortex-lattice  layout  which  overcomes 
these  objections  and  which  can  be  applied  to  a  wide  variety  of  configuration,  including 
rotary  wings. 
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RECENT  DEVELOPEMENTS  AND  SCOPE  OF  WORK 

(A)  2-D  thin  airfoil  theory: 

In  the  conventional  VLM.  a  thin  airfoil  is  divided  into  a  number  of  element  panels. 
N,  and  a  horseshoe  vortex  is  placed  at  the  1/4  quarter-chord  of  each  panel.  The  control 
point  is  located  at  the  3/4  quarter-chord  of  each  panel.  The  results  agree  with  James10 
analysis  that  the  first  prediction-value  (vortex  strength)  is  consistently  11.4%  too  low  for  a 
variety  of  cambers.  Furthermore,  increasing  the  number  of  elements  does  not  help  acurracy. 
But  C,  (lift  coeff.)  and  Cm  (pitching  moment  coeff.)  are  always  exact  for  N  greater  than  2. 
This  is  the  reason  that  the  majority  of  aerodynamists  use  conventional  VLM.  Lan3  developed 
the  so-called  "semi-circle  rule"  to  select  collocation  points  and  obtain  the  essentially  exact 
C|,  Cm  and  vortex  strenlh.  Kocuratc"  used  a  Doublet-Vortex  method  to  find  the  local 
circulation  in  this  2-D  case  and  applied  it  to  the  lifting-surface  performance  analysis  for 
hovering  rotors.  His  panel  spacing  is  biased  by  a  cosine  distribution,  but  control  points  are 
located  at  the  midpoint  of  each  panel.  The  method  does  converge  as  N  is  increased,  but  the 
rate  of  convergence  is  still  too  slow.  Results  in  this  paper  show  that  the  semi-circle  rule  can 
be  applied  to  both  panel  size  and  collocation  points  to  obtain  essentially  exact  values  of  Q, 
Cm,  and  vortex  strength. 

(B)  Lifting  line  theory: 

Glaucrti:  solved  Prandtl's  lif ting— line  formula  b\  using  a  Fourier  Series  method. 
Dejarnctli  applied  Multhopp's  intcpolation  technique  to  obtain  the  same  result  as  Glauerl. 


The  conventional  VLM  divides  the  entire  wing  span  into  a  number  of  panels.  M.  Also,  the 
trailing  point  is  located  at  the  boundary  of  each  panel;  and  the  the  control  point  is  placed  at 
the  midpoint  of  each  panel.  The  1/4-3/4  rule  is  also  used  by  Hough'4  to  afford  a 
singnificant  reduction  in  computational  costs.  Dejoung'5  has  proved  that  there  is 
mathematical  convergence  of  the  VLM  to  the  exact  answer,  when  the  1/4-3/ 4  rule  is  used 
as  M  (or  N,  in  2-D  case)  approaches  infinity.  The  crux  of  the  matter  lies  in  assessing  how 
fast  the  verified  1/4-3/4  rule  improves  the  convergence  and  if  it  is  an  optimized-choice. 
Actually,  the  1/4-3/4  rule  is  a  special  case  of  the  Finite-Difference  method,  for  which 
locations  presumably  can  be  chosen  arbitrarily.  In  this  paper,  however,  we  show  that  if  the 
collocation  points  are  laid  out  according  to  a  special  shape  function,  which  has  the  same 
shape  as  the  desired  unknow  circulations,  then  the  VLM  generates  the  fastest  convergence. 
These  shape  functions  must  have  quadradic  behavior  near  the  wing  tip  to  match  the 
asymptotic  tip  behavior16.  Both  parabolic  and  semi-circle  functions  are  valid  for  this 
requirement  at  the  tip.  and  these  functions  result  in  an  unequally -spaced  Finite-Difference 
technique.  Results  from  different  shape  functions,  compared  with  the  classic  solution  for  a 
elliptical  wing,  confirm  the  rapid  convergence  of  the  new  method.  A  mathematical  proof  (of 
the  optimality  of  the  correct  shape  function)  will  be  included  in  the  final  paper. 

Basin'7  derived  the  mathematic  equation  of  lifting-line  theory  for  a  rotary  wing  (in 
the  prescribed-wake  sense).  Rosen’  extended  it  to  curved  blades.  Both  of  these  use  10 
equally-spaced  meshes,  but  the  results  are  not  particularly  accurate.  Part  of  reason  stems 
from  the  fact  that  M  is  not  large  enough  due  to  poor  selection  of  collation  points.  By  using 
a  semi-circle  or  quater-circle  to  select  collation  points,  we  can  guarantee  rapid  convergence. 
Also  included  in  this  paper  is  an  extension  of  Mullhopp's  interpolation  technique  to 
Prandtl’s  lifting  formula  for  rotary-wing  problems.  This  technique  proves  to  be  equivalent 
in  accuracy  to  use  of  the  "Fourier  Series"  method.  Although  based  on  a  3-D,  rigid-wake 
concept,  an  example  is  given  for  lightly-loaded  hovering  blades,  calculated  by  both 
Finite-Difference  and  Fourier  Series  techniques.  These  are  compared  with  the  experimental 
data'",  and  both  give  excellent  correlation. 


(C)  Lifting-surface  theory: 

The  conventional  VLM  layout  for  lifting-surface  theory  is  to  use  uniformly-spaced 
panels  for  both  chordwise  and  spanwise  directions  over  the  whole  wing.  The  1/4-3/4  rule  is 
also  used  by  Hough14  for  the  spanwise  direction  and  by  Belotserkovskii6  for  the  chordwise 
direction.  A  reviews  of  these  results3-6  reveals  that  (M.N)  should  be  at  least  (40,2)  to  give 
good  convergence  for  typical  fixed  wings.  Dejarnette13  also  applied  Multhopp’s  interpolation 
technique  to  the  lifting-surface  formula  with  smaller  computation  time;  but  the  method  can 
be  used  only  for  rectangular  wings.  Some  criteria,  developed  from  combination  of 
lifting-line  and  thin— airf oil  theory,  can  be  applied  in  order  to  select  collocation  points. 
Results  in  this  paper  show  that  a  mesh  size  of  only  (10,2)  can  give  excellent  convergence  to 
exact  values  when  one  uses  the  semi-circle  rule  to  select  panels.  For  rotary  wings,  Rosen19 
and  Chang30  use  only  (10,1)  equally-spaced  panels  to  calculate  aerodynamics;  and  the  results 
seem  to  be  doubtful.  Proper  selection,  made  by  the  semi-circle  or  quarter-circle  rule,  should 
improve  the  convergence.  These  criteria  are  valid  not  only  for  VLM  but  also  for  the  Double 
Lattice  method,  which  is  an  extention  of  VLM.31  For  example,  in  the  case  of  the  propeller 
analysis  by  Murray33,  an  increase  in  outboard  unequally-spaced  panels  actually  appears  to 
have  a  detrimental  effect  on  convergence.  Thus,  meshes  cannot  be  selected  or  refined 
arbitrarily.  Finally,  one  can  reformulate  the  problem  to  solve  for  induced  drag  and  pitching 
moment  as  well  as  for  lift. 

CONCLUDING  REMARKS 

We  have  reviewed  the  conventional  VLM  literature  with  respect  to  selection  of 
collocation  points  and  the  corresponding  results.  These  generally  show  a  low  rate  of 
convergence.  A  new  criteria  for  choice  of  these  points  is  derived  from  lifting-line  theory.  A 
mathematical  proof  is  provided  to  show  that  this  is  the  optimum  choice.  In  summary,  the 
optimum  choice  of  collocation  points  can  be  done  as  follows: 

1.  The  choice  of  points  in  finite-difference  lift  calculations  should  be  made  on  the 


basis  of  equal  vortex  strength  between  points,  not  on  equal  spacing.  Ideally,  this  would 
require  an  adaptive  mesh  which  could  change  at  each  iteration;  but  practical  results  show 
this  is  not  necessary,  and  a  quasi-optimum  mesh  can  be  used. 


A 


2.  The  size  of  each  panel  must  smoothly  transition  from  large  increments  to  small 
increments,  moving  along  the  direction  in  which  the  rapid  change  of  vortex  stregth  occurs, 
and  following  a  smooth  shape  function.  This  function  must  be  quadratic  near  the  tip. 

3.  One  way  to  choose  the  points  in  a  quasi-optimum  manner  is  to  map  the  spanwise 
and  chordwise  interval  according  to  a  function  which  has  the  approximate  properities  of 
criterion  1  and  2,  such  as  a  semi-circle.  Both  panel  boundaries  and  collocation  points  must 
be  chosen  in  a  smooth  manner  according  to  this  function. 

4.  A  easy  way  to  map  the  interval  for  rotary-wings  in  subsonic  flow  is  to  use  the 
"Semi-circle  rule".  An  advantage  of  this  mapping  is  that  integrals  over  the  vortex  strength, 
can  often  be  reduced  to  closed  form.  Sevsral  semi-circles  can  be  used  if  several 
discontinuities  occur. 

Multhopp's  interpolation  technique  is  also  specified  as  a  semi-circle  rule  for 
unequally-spaced  Finite  Differences,  which  yields  the  same  error  as  the  Fourier  Series 
method  and  can  be  applied  for  rotary  wings.  The  results  for  wing  problems  by  the  Finite 
Difference  technique  (unequally-spaced  mesh)  are  found  to  compare  well  with  experimental 
data,  but  with  smaller  computational  times,  improved  accuracy,  and  simplfied  mathematical 
derivation,  as  compared  with  other  continuous  loading  methods. 
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(1)  FOURIER  SERIES  METHOD,  REF  12 

(2)  MULTHOPPS  INTERPOLATION. 
REF  13 

(3)  Panel  spacing  is  biased  by  a  parabolic 
distribution,  but  control  points  are  located 
at  the  midpoint  of  each  panel.  Also,  wing 
tips  are  specified  as  boundary  points. 

(4)  Panel  spacing  is  biased  by  a  coifsine 
distribution,  but  control  points  are  located 
at  the  midpoint  of  each  panel.  Also,  wing 
tips  are  specified  as  boundary  points. 

(5)  Both  control  and  trailing  points  are 
selected  according  to  equally-spaced 
panels  and  wing  tips  are  specified  as 
boundary  points. 

(6)  Panel  spacing  is  biased  by  a*  elliptic 
distribution,  but  control  points  are  located 
at  the  midpoint  of  each  panel.  Also,  wing 
tips  are  specified  as  boundary  points. 

(7)  Both  control  and  trailing  points  are 

selected  according  to  the  parabolic 

function  and  wing  tips  are  specified  as 
boundary  points. 

(8)  Both  control  2nd  trailing  points  are 
selected  according  to  semi-circle  rule  and 
wing  tips  are  specified  as  boundary  points. 

(9)  Both  control  and  trailing  points  are 

selected  according  to  equally-spaced 

panels  and  wing  tips  are  specified  as 
trailing  points. 

(10)  Both  control  and  trailing  points  are 

selected  according  to  the  parabolic 

function  and  wing  tips  are  specified  as 
trailing  points. 
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INTRODUCTION 


1.1  BACKGROUD 

Helicopter  vibration  reduction  is  important  as  well  as  useful.  With  less  vi¬ 
bration,  a  helicopter  can  survive  much  longer  and  can  provide  a  comfortable  en¬ 
vironment  for  passengers.  However,  unlike  conventional  fixed-wing  aircraft,  the 
helicopter  suffers  an  intrinsic,  severe  vibration  source  in  the  main  rotor,  which  con¬ 
tains  elastic  blades  and  is  connected  flexibly  to  the  fuselage  by  a  hub-pylon  system. 
So,  the  problem  is  rather  sophisticated.  It  is  well  known  that  the  fuselage  motions 
due  to  rotor  vibrations  will  cause  the  hub  to  move  in  all  degrees  of  freedom.  This 
hub  motion  can  cause  the  hub  loads  to  be  different  from  those  calculated  for  a 
fixed-hub  condition.  This  alteration  can  often  be  an  order-of-magnitude  change. 
Therefore,  in  studying  the  effect  of  hub  motions  on  hub  loads,  we  are  actually 
studying  a  feedback  system. 

1.2  PREVIOUS  WORK 

The  concept  of  performing  a  coupled  rotor/airframe  vibration  analysis  by 
impedance  matching  was  pointed  out  in  1964,  Reference  1,  which  indicates  two 
important  facts.  First,  a  coupled  rotor/ airframe  analysis  can  be  performed  in  a 
rigorous  maner  by  separate  calculation  of  rotor  and  fuselage  impedances  followed 
by  a  matching  of  forces  and  displacements  at  the  hub.  Second,  the  rotor  impedance 
need  only  be  calculated  for  a  single  blade  and  then  appropriately  transformed  to 
apply  to  any  number  of  blades. 

In  1974,  Staley  and  Sciarra  treated  the  vertical  vibrations  of  a  coupled  rotor 
and  fuselage,  including  the  effect  of  vertical  hub  motions.  They  used  a  lumped 
mass  for  rotor  impedance  and  showed  that  hub  motions  could  create  order-of- 
magnitude  changes  in  hub  loads.  In  Reference  3,  Hohencmser  and  Yin  further 
investigate  the  effects  of  rotor-body  coupling.  Their  model  for  rotor  impedance 
is  based  on  a  rotor  representation  tha.  includes  two  masses(each  equal  to  one- 
half  of  the  total  rotor  mass)  connected  by  a  spring  to  represent  the  first  flapping 


frequency.  Reference  3  presents  some  very  interesting  conclusions  that  pertain 
to  fuselage  design.  Particularly,  it  notes  that  under  certain  conditions  it  may  be 
desirable  to  tune  a  fuselage  frequency  to  the  blade  passage  frequency  in  ordor  to 
eliminate  hub  loads.  Also,  it  outlines  a  method  of  computing  the  complete  rotor 
impedance  by  finite  elements  and  transfer  matrices.  Other  work  on  the  impor¬ 
tance  of  hub  impedance  may  be  found  in  References  4-6. 

Since  rather  crude  models  have  been  used  for  hub  impedance  (rigid  mass,  no 
aerodynamics,  etc.),  one  might  wonder  why  more  sophisticated  models  were  not 
used.  The  answer  is  straightforward.  These  were  only  the  initial  investigations 
into  this  effect.  Furthermore,  although  most  analysts  realized  the  importance  of 
detailed  blade  modeling  for  fixed  hub  loads  (blade  modes,  unsteady  aerodynam¬ 
ics,  periodic  coefficients,  etc.),  it  was  not  clear  in  the  beginning  which  of  these 
effects  would  be  important  for  finding  the  role  of  hub  motion  on  loads.  Because 
of  the  high  frequencies  involed  (4/rev,  8/rev),  many  felt  that  inertial  terms  would 
dominate. 

Reference  7  offers  a  sophisticated  (but  linear)  rotor  flapping  model  that  allows 
for  a  detailed  investigation  of  both  loads  and  impedance  (even  in  the  presence 
of  periodic  coefficients).  The  method,  generalized  harmonic  balance,  involves  a 
computer-bassed  manipulation  of  equations  that  allows  many  degrees  of  freedom, 
many  modes,  and  many  harmonics.  In  Reference  8,  Hsu  and  Peters  apply  this 
method  to  a  flexible  rotor  and  then  use  impedance  matching  to  include  plunge, 
pitch,  and  roll  of  the  hub.  This  combined  technique  proves  to  be  very  efficient 
on  two  counts.  First,  the  calculation  for  only  one  blade  can  be  used  for  n-blades 
(as  in  Reference  1).  Second,  wholesale  changes  in  fuselage  properties  can  be  made 
without  a  requirement  to  recalculate  rotor  properties.  It  is  interesting  that  other 
investigators  who  began  with  a  full-blown,  coupled  analysis  later  changed  to  the 
impedance  matching  technique,  Reference  9-10. 

For  inplane  vibrations  with  a  fixed  hub,  Reference  8  sets  up  the  rotor  equa¬ 
tions  with  rigid  blades,  and  Reference  12  treats  the  same  problem  but  with  elastic 
blades. 


In  Reference  13,  a  fuselage  model  with  offsets  and  9  degrees  of  freedom  is 
considered,  coupled  with  a  rigid  inplane  blade  model.  There  is  an  important  con¬ 
clusion  in  reference  13  which  will  be  very  useful  for  further  nonlinear  analysis. 
It  is  that  the  addition  of  inplane  degrees  of  freedom  does  not  significantly  affect 
the  plunge  vibrations  for  the  cases  considered,  and  these  cases  are  for  resonable 
configurations. 

1.3  SCOPE  OF  WORK 

In  order  to  predicate  the  helicopter  coupled  vibration  much  better,  the  rotor 
with  elastic  blades  should  be  considered.  The  scope  of  this  research  is  to  continue 
the  previous  efforts  in  the  study  of  inplane  coupled  rotor-body  vibrations.  To 
do  so,  first,  the  dynamic  partial  differential  equations  of  elastic  blades  including 
hub  motions  should  be  derived.  Second,  hub  loads  equations  should  be  derived. 
Third,  the  fuselage  equations  are  taken  to  be  the  same  as  those  in  Reference  13. 
The  major  goal  of  this  research  is  to  investigate  if  higher  frequencies  of  the  elastic 
blades  will  affect  the  coupled  rotor-body  vibrations. 


-6- 


2.  ELASTIC  BLADE  EQUATIONS 


The  equations  of  motion  for  a  flexible  hingeless  rotor  blade  in  forward  flight, 
as  derived  in  Reference  12,  are  the  starting  point  for  this  analysis.  The  present 
analysis,  however,  expands  the  previous  work  in  following  way.  The  hub  motions 
(including  plunge,  longitudinal  and  lateral  displacements,  pitch  and  roll  angles, 
and  their  first  and  second  derivatives)  are  added  to  the  equations. 

2.1  COORDINATES  AND  THEIR  TRASFORMATIONS 

The  coordinate  systems  used  in  the  present  analysis  are  shown  in  Figure  1. 
The  triplet  X,Y,Z  represents  a  fuselage  fixed  coordinate  system;  and  the  triplet 
x,y,z  represents  a  rotating  coordinate  system.  Between  the  Z  and  z  axes  there  is 
a  precone  angle.  The  deflection  components  of  the  elastic  axis  of  the  deformed 
blade  (u,  v,  and  w)  are  taken  in  the  x,y,z  coordinate  system.  The  final  set  of  axes, 
x‘,rS,z'  are  taken  along  the  deformed  axes. 

There  are  transformations  among  these  coordinates.  Figure  2  shows  the  pitch 
angle  ac  and  roll  angle  a,;  and,  Figure  3  indicates  azimuth  angle  4>  and  precone 
angle  0pc.  With  these  angles,  one  can  trasform  any  vector,  say  F,  from  the  fixed 
coordinate  X,Y,Z  to  the  undeformed  rotating  coordinate  x,y,z.  We  have: 
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The  transformation  from  undeformed  coordinates  to  deformed  coordinates  is  given 

by 


where, 


(3) 


(  1 

[Td\  =  -9w+  -  v+ 

^  —  tu+  +  6v+ 


—9  —  w+v+ 


(^) 


With  these  transformations,  we  can  transfer  the  vectors  to  any  coordinates  needed. 
Other,  kinematic  contributions  to  9  are  treated  as  pitch-lag  coupling. 


2.2  HAMILTON’S  LAW 


The  next  step  in  the  derivation  is  to  apply  Hamilton’s  Law  to  obtain  equations 
of  motion.  The  Hamilton’s  Law  can  be  expressed  as 


6T~6W)dt+^6qi 

aqi 


=  0 


(5) 


where,  SU  and  ST  are  the  variations  in  strain  energy  and  kinetic  energy,  and  is 
the  virtual  work,  ‘g,-’  is  the  generalized  coordinate.  To  apply  equation  (5),  SU,  ST 
and  SW  must  be  formulated  in  terms  of  the  generalized  coordinates. 


2.2.1  STRAIN  ENERGY 


The  strain  energy  can  be  expressed  as  follows: 


U  =  r  j  JJ  Ei^xdA  dx 


then, 


SU 


i:il 


EexxScxx  dA  dx 


(6) 


(7) 


v"\t)  cos  9t  —  f  cos  9e\  —  w"[?7  sin  9e  +  j  cos  9e] 


(«) 


where, 


Therefore,  SU  can  be  fuether  expressed  as 


SU  =  J  {vx.\6u'  +  v'6v'  +  w'Sw']  +  [Mz. +  M„>Re6]6V 
+  [MX.R<6  -My.]6W“}dx 


where, 


f  ,  v’^  tu,a  1 

VX,  =  EA{«'  +  T+-T'j 

Mgl  =  EIy‘[v"Red  -  w"\ 


Mg'  =  EI,\v"  +  w"RJ] 


Integration  by  parts  of  SU  gives  the  following  result: 


rR 

6U~  (?u6u  +  ?v6v  +  ?v,6w)dx  +  b(U) 

Jo 


where, 


K  =  -(Vg,y 


?v  =  \M,>  +My'Re9\"  -  (Vgtv')' 
ya  =  [M,>Re9  -  My']"  -  (Vg'w'Y 


and  the  boundary  condition  is 


6({/)  =  Kx«5u|  +  jv,V  -  [M,.  +My.R'd]'  j^| 

+  [M,<  +My>Re6\5v'\  +  { Vx.w  -  [Mg.Rtd 

lo  1 

\* 

+  [Mg.Re9  -  My.]Sw'\ 


(18) 


2.2.2.  KINETIC  ENERGY 


The  kinetic  energy  can  be  expressed  as 


- UIL  pV  -V  dr}  d{  dx 


and  the  variation  is 


- ai 


pV  ■  SV  drj  df  dx 


where  the  velocity  vector  is 


V  =  tjz’i  +  dazi  sin  ip  +  acZ\  cosip  -  a, yiftpC  cos  ip 

+  acyif)pcsmip  -  Qyi  -  X  cos  ip  +  Y  sin  ip  -  Zfipc  J 
+  /{jM  +  dt2i  cos  ip  -  aczi  si nip  —  Clzifdpc  +  X  sin  ip 
+  d,x\0pc  cos  ip  -  dcXiffpc  sin  ip  +  0^!  +  Y  cos^j 
+  £{*i  -  d,y !  cos  ip  +  acy\  sin  ip  +  Uyi0pc  —  a,xi  sin  ip 
—  dex  1  cos  ip  +  X0pc  cosip  —Y  fipc  sin  ip  -  Z  J 


and  where, 


Xi  —  x  +  u  —  v'\r)  cos  9  —  f  sin  0]  —  w'\q  sin  8  +  f  cos  6) 


y\  =  w  +  rj  cos  9  —  f  sin  9 


ij  =  to  +  rj  sin  9  +  f  cos  9 


After  integration  over  the  blade  cross  section,  the  variation  of  kinetic  energy  be¬ 


comes: 


r « 

ST  =  I  {ZM  +  Zv6v  +  Z„6w  +  Zv,6v'  +  2u.6rf)  dx 

Jo 


where, 


Zu  =  m  [«  —  dgiii  sin  fit  +  2ni  +  X  cos  Of  —  Y  sin  f)£  —  D acw  sin  f)< 


+  2Qdg0pcx  cos  Ctt  -  2 Qdc0pcx  sin  Qt  +  Cldaw  cos  Ht 


ojf*  fmfmjCmjCmjtmjtrnJl ZmJl 


>•*.  VN*  »*  ’  *V  .*  *.*  'a  Vs.'  *  J 


+  fi2a;  —  a,ti/sinOf  +  a',2s  +  2  dtdcx  cos  fit  sin  fit 


—  acw  cos  fit  —  Z0pc  +  Zat  sin  fit  4-  Zac  cos  fit] 

+  m( w  +  e  sin  0)[— a,  sin  fit  —  fit*',  cos  fit  —  Vli9pc\ 

Zv  =  m  [— v  +  2 dcw  sin  fit  +  0^pcv>  —  d,9pex  cos  fit  +  o', U0pcx  sin  fit 

+  QcPpcX  sin  fit  -+  dc9pcx  cos  fit  —  X  sin  fit  —  Y  cos  fit  —  2 d,w  cos  fit 
+  o’,2zsinfit  cosflt  —  acalp'haa(cos3nt  —  sin2 fit)  —  o'c2xsinfit  cos  fit 
+  Cl0pcw  —  Qde0pcx  sin  fit  —  ndcftpcx  cos  fit  +  Zd,  cos  fit  —  Zdc  sin  fit] 
+  mfi2  (i>  +  ecos  8)  —  2Cimu  +  2mefi(v'  cos  8  +  w'  sin  9) 

+  m(w  +  esin0)[— a,  cos  fit  +  dc  sin  fit] 

Zu  =  m  [— w  +  2 d,v  cos  fit  —  2 dcv  sin  fit  —  2fi 9pcv  +  a,x  sin  fit  +  2a‘,fix  cos  fit 
+  acx  cos  fit  —  2fi acx  sin  fit  +  Y  d„  +  X0pc  cos  fit  —  Y  9 pc  sin  fit 
—  Z  —  Xalphac  —  CPfipcX  ~  2ClX0pc  sin  fit  —  2CIY 9Pc  cos  fit 
+  m(v  +  e  cos  0)  [a,  cos  fit  —  2fia',  sin  fit  —  dc  sin  fit  —  2fio'f  cos  fit] 

Zv'  —  me  cos  0[Y  sin  fit  —  X cos  fit  —  fi2x  —  2fWj 
z„.  =  me  sin  8[Y  sin  fit  —  X  cos  fit  —  fi2x  —  20i] 


in  which, 


m  =  JJ  p  dr]  d< 

(29) 

me  cos  9=  JJ  p(yi  —  v)  dr/ df 

(30) 

me  sin  9  =  JJ  p(zx  —  w)  dr)  df 

(31) 

Integration  by  parts  again  gives 

f  Zv>8v'dx  =  f  2vi6(6v) 

Jo  Jo 

ifl  rR 

=  Zv'8v\  -  Z'v,6v( 

n  JO 


and 


■N. 
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fR  lft  fR 

I  2„>Sw'  dx  =  2a'Sw\  —  I  2ui6wdx 
J  o  lo  Jo 

Finally,  the  variation  of  kinetic  energy  can  be  expressed  as 

rR 

6T  =  I  [2u6u  +  ( barZv  —  2'v,)6v  +  [2a  —  2'ul,)Sw ]  dx  +  b(T ) 
J  o 


(33) 


(34) 


where  the  boundary  condition  is 


b{T )  =  2v.Sv 


+  Z„iSw 


(35) 


2.2.3  VIRTUAL  WORK 


The  final  step  necessary  to  compute  the  equations  is  the  computation  of  virtual 
work.  The  virtual  work,  8W ,  is  mainly  due  to  the  nonconservative  external  forces 
which  come  from  the  theory  of  aerodynamics.  The  virtual  work  can  be  expressed  as 


8W  —  f  [L„6u  +  L„Sv  +  LuSw)  dx 

Jo 

where  Lu,  Lv,  La  are  forces  acting  along  x,y,z  directions. 


2.2.4  GENERAL  EQUATIONS 


(36) 


The  functional  of  Hamilton’s  Law  has  three  variables  «,  v,  and  w.  This  will 
lead  to  three  Euler-Lagrange  equations  which  are  the  required  equations  of  mo¬ 
tion.  Later,  u  can  be  eliminated  which  will  reduce  the  system  of  equations  to  only 
two.  The  basic  equations  come  from  Equation  (5).  If  we  ignore  the  trailing  terms 
that  are  cancelled  by 


then  we  have 


(37) 


Therefore, 


6T-6W)dt 


=  0 


(38) 


SU-ST-SW  =  0 


(39) 


Vy  .• 
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with  trailing  terms  dropped.  In  the  case  considered,  we  have 

?R 

J  |(f0  -  2tt  -  Lu)Su  +  [?„  -  (2„  -  2's)  -  Lv]Sv  4  [?»  -  {2a  -  2‘m.)  -  Lm]  Sw  }  dx  +  6(17)  -  b(T) 

-0  (40) 

From  the  principles  of  variational  calculus,  the  equations  of  motion  and  boundary 


■ 

conditions  can  be  obtained  from  Equation  (40). 

SU  equations: 

. 

Fy  —  2xi  -  6y 

(41) 

6V  equations: 

Yv  —  (2V  —  2vi)  =  L„ 

(42) 

•  , 

SW  equations: 

tv 

Y„  —  ( 2U  —  2ai )  =  La 

(43) 

and  boundary  condition: 

b{U)  -  6(T)  =  0 

(44) 

2.2.5  BASIC  EQUATIONS 


With  the  help  of  the  ordering  scheme  in  Reference  12,  but  keeping  flapping 
equations  o(e),  lead-lag  equation  o(«-3),  and  eliminating  axial  deflection  from  SU 
equation,  we  obtain  flapping  and  lead-lag  equations  in  nondimensional  form  as 
follows: 

Flapping  equation: 


A?te+  +++  +  (Aj  -  A3).R«0S+++  +  -  r  +  w  +  -  rw++  4»"  -  2v*+  -  c,*t*i  sin  p 
—  2o£icos  ip  —  orj*  cos  ip  4  2a* £  sin^  4  X**0pccos  ip  —  Y**  0pc  sin  ip  —  2**  4  9pci 

=  La  (45) 

Lead-lag  equation 

Ajs++  +  +  4  (Aj  -  A?)i?,06aru;++++  -  t+v+  -  tv++  -  v  -  2  /  w  +  w*+  dx 

J  o 

4  H**  —  a**  0pc  sin  ip  —  a**xw  sin^»  4  a*J£  sin  ip  cos  ip  —  2a*  w *  sin  ip  4  22*  a*  sin  y 
4  a\*20pc  cosip  4  aj*w  cos  ip  —  a**£  sin  ip  cos  ip  4  2a*tu>*  cos  ip  —  22*  a*t  cos  ip 
4  aJa*2(cos3^  —  sin3^)  —  20pcw*  4  X**  sin  ip  4  Y**  cos  ip 
=  Lv  (46) 


_  'w 


->  *  *.  _S  . 


2.3  AERODYNAMIC  FORCES 

In  deriving  the  aerodynamic  forces,  the  following  assumptions  are  made 

1.  Linear,  quasi-steady  aerodynamics 

2.  Blade  stall,  compressibility  and  reversed  flow  are  neglected. 

3.  Uniform  induced  flow 

Figure  4  shows  the  free  airflow  passing  through  the  helicopter.  And,  Figure  5 
shows  the  blade  element  geometry.  The  total  velocity  with  respect  to  air  mass  in 
forward  flight  will  be 


UT  =  Uair  -  U 


where, 


0air  =  Voof  -  (K  +  v)K  (48) 

and  V  is  the  velocity  obtained  in  Equation  (21).  Using  the  transformation  [Tjt] 
and  \Ti\  in  Equations  (2)  and  (4),  We  obtain  the  required  velocity  components  in 
deformed  coordinate  with  dimensionless  form  as  follows: 

Dy‘  **  —Otii+ficos  ip  —  ii+  fi  cos  ip  —  n  sinip  —  v*  —  barx  —  u 

—  cos^  +  a^id  sinip  +  iD0Pc  —  ct*20pc  cos  ip 

+  sin  ip  —  X *  sinip  —  P*  cos  ip  —  9 A 

—  90pc(i  cos  ip  +  9aen  —  Bui  *  +  Oa,*x  sin  ip 

+  9a*i  cos  ip  +  02*  —  \ac  sin  ip  +  \a,  cos  ip  (49) 

€,>  =  —w+/icos  ip  —  w+X *  cos  ip  +  iD+f*  sin  ip  —  A 
+  Opsin  ip  +  9i  +  X*6  sin  ip  +  Y*6  cosip 

—  nf3pc  cos  ip  +  nae  -  w*  -)-  a't2  sin  ip 

+  a*£cos  ip  —  X*0pc  cosip  +  Y*  f)pcs\n  ip  +  2*  (50) 

In  the  deformed  coordinate,  the  dimensionless  aerodynamic  forces  can  be  expressed 


F,>  =  0 

F,'  =  j[«.  - 

f,,  =  -  J|0,.tV) 


•  *  «  "  «  **  .  '  *  “  -  ’  .  1  «,  •  .  «  .  <.  « 


V  vv  ».*  •  ’.C 


✓  V.  / 


Then,  we  can  transfer  the  aerodynamic  forces  above  to  undeformed  rotating  co¬ 
ordinates,  which  are 


Then,  the  external  aerodynamic  forces  are 


Lv=Fy 


2.4  FINAL  EQUATIONS 

Combining  Equation  (45),  (46)  and  Equation  (55),  (56),  we  obtain  the  com¬ 
plete  flapping  and  lead-lag  equations  as  follows: 

FLAPPING  EQUATION 

AjtD++++  +  (A3  -  A.\)R'8v  ++++  -T+w+  -  TW++  +  a**  -  2V*+  -  a*,*x sin  tp 
—  2a\x  cos  ip  —  a**xcos^  +  2a*2sin^  +  X**0pccos  tp  —  Y**0pc  sin  Tp  —  2 **  +  0pcx 

=  ^  |^/i3sin2  tp  +  6x'1  +  20pxs'mtp  +  M3acsin^  +  /«**£  sin3  tp  +  f*a*X  sintpcostp 
+  1*2*  smtp  +  i*aci  —  w*X  +  a£23  sin  tp  +  a*£3  cos  tp  —  u/+n3sin  tp  cos  tp  +  z*  i 
—  n  Asin  tp  —  \x  —  u7Ppc  sin  tpcos  tp  —  ftw*  sin  tp  —  w+x  cos  tp  —  n0pcX  cos  tp  j  (57) 


LEAD-LAG  EQUATION 


A3e++++  +  (A3  -  A?)i?e0tB++++  -  r+v+  -  r f?++  -  v 
—  2  f  w+w*+  dx  -f-  v**  —  a**0pC  sin  tp  —  a**2w  sin  tp 

J  o 

+  QfJ2isinV'cosV’  —  2or*tS*  sin  tp  |-  2  Z*a*  sin  tp  +  a*a*x0pccos  tp 
+  Ot*t*tjj  cos  tp  —  ot**x  sin  tp  cos  tp  -f-  2ajti>*  cos  tp  —  22* a*  cos  tp 


+  a*arj2(cos3^  —  sin3^)  —  20pcw*  +  X**  sin  tp  +  P**  cos  tp 


|tC+a/iacos2V>  +  A2  +  fi1 0,e7  cos7  ip  +  ft7a7e 

+  tE*2  +  Z*7  +  Qj323sin2^  +  a*7  £7cos7ip 

—  ib+9fi7  cos  ip  sin  ip  —  iD+ 9  ft£  cos  ip  +  2XiI>+ fi  cosip 
+  2tZ>+ fi7  f)pe  cos  ip  —  2w+ n7ae  cos  ip  +  2w+ w*  fi  cos  ip 

—  2iD+ a*tfi£  cos  ip  simp  —  2&+ a* picos7  ip  —  2w+  2*  ft  cos  ip 

—  9  Xfi  sin  ip  —  9 n7  ft pc  sin  ip  cos  'p  +  9ft7ac  sin^ 

—  H)*9fisinip  +  9a\nxs\r?ip  +  9a* fix  simp  cosip 

+  92* ft  simp  —9 Xi  —  9fif}pcx  cosip  +9fiacx  —  9w*x 
+  9 Z*£  +  9a*x 7  sin  ip  +9  a*  £2  cos  ip  +  2A nflpc  cos  ip 

—  2A fiac  +  2AiZ>*  —  2Aaj2  sin  ip  —  2Aq*2  cos  ip  —  2X2* 

—  2fi7fipcac  cos  ip  +  2iB*fiPpccos  ip  —  2a*j}pcfi£sinip  cos  ip 

—  2a*0pcfi£cos7ip  —  2s*  fi(ipe  cosip  —  2w*  fiac  +  2arca*2sin  ip 
+  2a* a,fi£ cos  ip  +  2 2*(iac  —  2 w*a*t£  sin  ip  —  2 w*a*£  cos  ip 

—  22*  ib*  +  2a£a*2a  sin  V1  cos  ^  +  22*  a\£simp  +2Z*a*2cos^> 

- ^[2a  +  p7sin7  +  2 fix  sin  ip  +  2fix*sin7ip 

+  2xX*  sin  ip  +  2fi?*  sin  ip  cos  ip  +  2iY*  cos  ip 

THE  BOUNDARY  CONDITION 


Vx>Su 


+  {Vx,v’  -\M,.+ My’ R<ey}6v 


+  [M'.+MylR'9}8  vr 


+  {Vx.w  -  \M,.Rt9  -  My.]'}6w  +  \M,'Rt6  -MyjW 


—  me  cos  9[Y  sin  (It  —  X  cos  (It  —  O7 x  —  2fl v\6v 


—  me  sin  9[Y  sin  (it  —  X  cos  (it  —  (i3x  —  2Qv\6w 


>0 


(58) 


(59) 


It  may  be  noted  that  for  a  rigid  blade  (for  which  w  =  px  and  v  =  f x ,  the  flapping 
Equation  (57)  and  the  lead-lag  Equation  (58)  reduce  to  those  of  Reference  13  with 
the  exception  of  nonlinear  higher-order  terms. 


3,  HUB  LOAD  EQUATIONS 

The  hub  loads  including  the  shears  and  moments,  in  both  rotating  and  fixed 
coordinate  systems,  are  shown  in  Figure  6.  Cx ,  Cy,  and  C,  are  radial,  lead-lag, 
and  flap  shear  components;  and  C„,  CB,  Ca  are  torsional,  in-plane  and  out-of-plane 
moment  components  in  the  rotating  coordinate  system.  Cz ,  Cx ,  and  C\  are  the 
vertical  thrust,  forward  drag,  and  side  force  components  in  the  fixed  coordinate 
system.  The  loads  with  respect  to  the  rotating  system  are  called  rotating  loads 
and  those  corresponding  to  the  fixed  system  are  called  nonrotating  loads.  They 
are  derived  in  detail  below. 

3.1  ROTATING  LOADS 

The  expressions  for  the  root  shear  and  moment  components  in  the  rotating 
coordinate  system  are  formed  by  spanwise  integration  of  the  aerodynamic  and  in¬ 
ertia  loads  as  follows: 


CO  1  f~ 

1! 

C  FxdX~-  [\oxdi 

Jo  7  Jo 

(60) 

c.-i 

7. 

fo  Pyd£~li  a°ydi 

(61) 

C,=  - 
7 . 

f  F,  dt  —  —  f  a0,  dx 

Jo  7  Jo 

(62) 

Cut  —  ~ 
7 . 

fl  3  cl  3  y i  3  r1 

/  F,xdx  +  -I  a0,idx  +  -l  Fxwdi - /  aoxwdx 

Jo  7  Jo  7  Jo  7  Jo 

(63) 

Cu  «0 

(64) 

C„  S3  0 

(65) 

where  the  aerodynamic  forces  are  obtained  in  Equation  (54),  and  the  dimention¬ 
less  accelerations  are  as  follows: 


=  a**  -  2a* w*  sin  ip  +  2a*cw*  cos  ip  -2v*  -  a*Jisin2^ 
-  a*7t  cos  ip  +  2a*a*t  sin  ip  cos  ip  -  2a*a0pcx  cos  ip 
+  2 a*0pc£  sin  ip  -  2  -  w  4-  dpc®  ~  a**w  cos  ip 


+  al* iD  cos  $  + X**  cos  $+¥**  sin  i>  (66) 

Boy  =  V**  +  2fl*  —  2a* t fl*  cost/)  —  2a*  ID*  sin^  —  20pcil* 

—  a*tw  sin  ip  +  aju)  cos  ip  —  «  —  aj32  sin  ip  cosip 

—  a*a*t& sin3  ip  +  aja*2  cos3  ip  +  a\$pc5,  sin  ip 

+  a*3i  sin  ip  cos  ip  +  aypipci  cos  ip  -  a**0pcx  sin  ip 

—  a**ti>  cos^  +  ajtesin  ip  —  aj*u>  sin^>  —  a*wcos  ip 

—  X**  snnip  +  P**  cos  ip  (67) 

B0,  —  »**  f  2a*  2  cos  ip  +  2a*  £  sin  ip  +  ppcZ 

+  a**±  sin  ip  -  aj*s  cos  ip  +  2**  (68) 

3.2  NON-ROTATING  LOADS 

The  shear  and  moment  components  in  the  fixed  coordinate  system  are  formed 
from  the  load  components  in  the  rotating  coordinate  system.  They  can  be  simply 
expressed  as 

Cx  —  C,  cos  ip  —  Cs  sin  ip  (69) 

Cy  —  Cj  sin  ip  +  Cy  cos  ip  (70) 

Cz  =  C,  (71) 

Cl  —  ~C„  sin  ip  (72) 

Cfu  =  Cv  cos  ip  (73) 

They  can  be  put  into  matrix  form  as  done  in  Reference  13.  It  is: 


4.FUSELAGE  EQUATIONS 


In  this  research,  the  same  model  as  in  Reference  13  is  used  for  fuselage,  except 
that  X  and  Y  axes  have  opposite  positive  directions.  Figure  7  shows  the  fuselage 
model  in  the  longitudinal  and  pitch  dii  actions. 

This  model  includes  9  degrees  of  freedom.  These  are: 

1)  vertical  rigid-body; 

2)  rigid-body  pitch; 

3)  rigid-body  roll; 

4)  rigid  lateral; 

5)  rigid  longitudinal; 

6)  elastic  vertical; 

7)  elastic  lateral; 

8)  elastic  pylon  in  pitch; 

9)  elastic  pylon  in  roll. 

The  model  also  includes  vertical  offsets  between  the  fuselage  center  of  mass,  the 
pylon  focus,  the  pylon  center  of  mass,  and  the  rotor  center. 

The  fuselage  is  modeled  as  a  uniform  beam  with  a  lumped  mass  A/r  added 
at  the  center.  The  mass  of  the  pylon  is  separated  from  the  fuselage,  which  is 
connected  to  the  fuselage  through  pitch  and  roll  torsional  springs.  The  fuselage 
equations  of  motion  are  the  same  as  those  in  Reference  13,  except  that  the  X  and 


Y  axes  are  reversed. 


r 


5.  SOLUTION  METHOD 

5.1  GALERKIN  METHOD 

The  flap-lag  equations  of  motion,  Equation  (57)  and  (58),  are  nonlinear,  vari¬ 
able  coefficient,  integro-partial  differential  equations.  First,  these  equations  can 
be  reduced  to  nonlinear  ordinary  differential  equations  in  terms  of  the  generalized 
coordinates  by  use  of  the  modal  form  w  =  and  v  =  ^p(,  and  the  Galerkin 
method  as  below: 


Where  M,  <7,  and  K  are  functions  of  q ,  p,  and  their  derivatives  and  where 
*,/  =  l,...m  number  of  flap  degrees  of  freedom; 

and 

i,J  =  l,...n  number  of  lead-lag  degrees  of  freedom. 

5.2  HARMONIC  BALANCE 

As  done  in  Reference  13,  fuselage  equations  can  be  transformed  into  a  set  of 
linear,  algebraic  equations  by  a  harmonic  balance  which  is  formulated  in  matrix 
notation  for  a  linear,  multi-degree-of-freedom  system  in  Reference  7.  In  our  case, 
the  equations  are  nonlinear.  Therefore,  an  important  conclusion  obtained  in  Ref¬ 
erence  13  can  be  used.  As  mentioned  before,  it  says  that  the  addition  of  inplane 
degrees  of  freedom  does  not  significantly  affect  the  plunge  vibrations  for  the  cases 
considered,  and  these  cases  are  for  resonable  configurations.  Figure  8  and  9  shows 
the  comparison  of  plunge  vibration  with  or  without  inplane  degrees  of  freedom. 

Therefore,  we  can  first  solve  the  plunge-pitch-roll  problem  with  inplane  pa¬ 
rameters  eliminated.  To  do  this,  the  flapping  equation  is  expressed  as  below: 
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and  the  load  equation  (74)  becomes 


-[*m]  +[*{*)]  «.  +[bw]  «;  +  [%-)]  a;* 
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For  the  fuselage  equations,  we  have  eliminated  X,  F,  and  ?F,  and  their  derivatives. 
Now,  these  equations  are  linear  ordinary  differential  equations  with  peoriodic  co¬ 
efficients.  They  are  easily  solved  by  harmonic  balance  technique  and  impedance 


matching.  The  control  vector  components  fft,  6e,  A,  /?,  0pc,  and  ac  can  be 
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obtained  from  a  trim  procedure. 

For  the  next  step,  we  go  to  Equation  (75)  with  gi  ...qj  and  acy  a,,  and  Z  and 
their  derivitives  as  knowns.  The  flap-lag  Equations  (75)  can  be  further  expressed 

as 


and  the  load  Equations  (74)  can  be  further  expressed  as 


Then,  using  flap-lag  Equations  (78),  load  Equations  (79)  and  complete  fuselage 
equations,  we  can  finally  solve  the  inplane  problem. 
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6.  STATEMENT  OF  WORK 
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In  helicopter  preliminary  design  stage,  it  is  very  important  for  designers  to 
determine  the  helicopter  structure  dynamic  parameters.  Besides  getting  a  stable 
system,  we  hope  that  helicopter  vibrations  can  be  reduced  as  much  as  possible 
in  the  helicopter  flight  envelope.  Since  different  combinations  of  fuselage  and  ro¬ 
tor  parameters  can  cause  quite  different  vibration  levels,  we  should  study  them 
carefully.  As  is  well  known,  most  hingeless  rotors  can  be  considered  as  relatively 
soft  flapwise,  except  the  ABC  blade  which  is  stiff  flapwise.  So,  we  will  choose 
a  soft  flapwise  rotor  system  as  done  in  Reference  7  and  Reference  13.  For  the 
inplane  direction,  soft  inplane  and  stiff  inplane  both  are  common.  Therefore,  this 
investigation  includes  both  soft  inplane  and  stiff  inplane  rotors.  In  brief,  we  will 
change  fuselage  dynamic  parameters  (the  fuselage  first  bending  natural  frequen¬ 
cies  in  vertical  and  lateral  directions,  and  the  pylon  torsional  natural  frequencies 
in  pitch  and  roll  directions)  to  match  two  categories  of  rotor  systems,  which  are 
soft  flapwise,  and  both  soft  inplane  and  stiff  inplane.  Also,  the  fuselage  layout 
will  affect  the  dynamic  response  significantly,  as  obtained  in  Reference  13,  the 
configurations  without  or  with  offsets  are  concerned. 

In  order  to  finally  solve  the  Equations  (78)  and  (79)  with  different  dynamic 
parameters  described  above,  the  detailed  procedure  are  as  follows.  At  first,  Equa¬ 
tions  (76)  and  (77)  should  be  solved  in  order  to  obtain  the  generalized  coordinates 
9i  ...qj  and  ac,  at,  and  2  and  their  first  and  second  derivatives  as  knowns.  To  do 
so,  three  nonrotating  out-of-plane  bending  modes  are  used,  which  are  the  same 
as  those  in  Reference  8.  Also,  the  baseline  parameters  are  chosen  as  the  same 
as  those  in  Reference  8  for  comparison  purposes.  The  baseline  parameters  are  as 
follows: 


Rotor: 

4  blades 

1  =  8 

/*  =  •  3 

p=  1.12 

wj  =  2.5 

W3  =  4.5 

Fuselage 

f gm  =  -37 

fgL  =  -14 

Qfv  =  1.53 

Ufm  =  2.58 

W/t  =  l-lSi'/m 

O/v/wew  =  1-45 

&Jm/ucm  =  10.0 

ujlIucl  —  4.47 

C4 
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II 
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to 
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tSS 

II 
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to 

O 
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II 

II 

II 

Trim  Condition: 

\CZ\  =  .0144 

\Ci\  =  \Cm\  =  0.0 

Oc  =  -.0715 

\  =  .03187 

•  •  ^  . 


Next,  we  add  three  nonrotating  blade  inplane  bending  modes.  In  order  to 
study  the  dynamic  response  for  both  soft-inplane  and  stiff-inplane  rotors,  two  sets 
of  three  nonrotating  blade  inplane  bending  modes  are  used  to  correspond  to  dif¬ 
ferent  rotors.  For  the  comparison  purposes,  the  first  inplane  nonrotating  bending 
frequencies  corresponding  to  soft  inplane  and  stiff  inplane  rotors  are  the  same  as 
those  in  Reference  13,  which  are  wf  =  .7  and  ut  =  1.4,  respectively.  Also,  we  need 
both  lateral  damping  coeffecients  of  the  fuselage  and  the  steady  portion  of  side 
forces,  which  are  gy  =  .02,  .002,  and  1(7x1  =  |Cv|  =  0.0,  the  same  as  those  in  Reference 
13. 

The  last  step  is  to  study  the  effects  of  offsets.  The  same  parameters,  h  =  .4, 
ir  =  0.0,  and  if  —  .2  as  those  in  Reference  13,  are  used  for  both  soft  inplane  and  stiff 
inplane  cases.  Thus,  we  can  obtain  all  responses  for  different  possible  situations. 

At  the  present  point  in  our  research  we  are  transfering  our  computer  codes 
from  the  VAX  system  at  Washington  University  to  the  CDC  system  at  Georgia 
Institute  of  Technology.  When  this  work  is  completed,  we  will  be  able  to  continue 
the  numerical  work  and  should  have  our  first  answer  shortly.  The  first  runs  will 
be  validation  runs  to  ensure  that  we  can  repeat  the  results  of  References  8,  12, 
and  13.  Then  we  will  procede  to  studies  of  parameter  variations. 


7.  NOMENCLATURE 


a 

=  slope  of  lift  curve,  l/rad 

A 

=  ratio  of  rotor  mass  to  moment  of  inertia,  MxR/lv 

b 

=  number  of  blades 

Of. 

=  blade  profile  drag  coefficient 

6z 

=  conventional  thrust  coefficient,  thrust /pv&R4 

Cx,Cy,Cz,Cm,Ci 

=  vibratory  portion  of  nondimensional  longitudinal 

force,  lateral  force,  thrust,  pitch  and  roll 

moment  over  ca 

Cx,Cy,Cz,Cm,Cl 

—  steady  portion  of  coefficients 

if 

=  offset  between  focus  and  center  of  pylon, 

divided  by  R 

*p 

—  offset  between  focus  and  center  of  pylon, 

divided  by  R 

ir 

=  offset  between  hub  and  center  of  pylon, 

divided  by  R 

t 

=  mass  centroid  offset  from  elastic  axis  of  a  blade 

El 

=  beam  cross-section  bending  stiffness 

{?) 

=  vector  of  forces 

0 

=  nondimensional  acceleration  of  gravity,  g/Q^R 

QtiQyiQmiQl 

=  plunge,  lateral,  pitch  and  roll  structural 

damping,  «  2f/w„ 

h 

=  offset  between  hub  and  focus,  divided  by  R 

im 

=  fuselage  receptance 

w 

=  identity  matrix 

ipt 

=  pitch  inertia  moment  of  pylon,  divided  by  A fpR'1 

=  roll  inertia  moment  of  pylon,  divided  by  MPR2 

Ar» 

=  pitch  inertia  moment  of  fuselage,  divided  by  Mr  R 

I*, 

=  roll  inertia  moment  of  fuselage,  divided  by  A/f/?3 

l 

=  length  of  the  beam,  m 
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=  mass  per  unit  beam  length,  kg/m 
=  lumped  mass  on  the  center  of  fuselage,  kg 
=  mass  of  pylon,  kg 
=  mass  of  fuselage,  Mc  +  ml,  kg/m 
=  mass  of  whole  fuselage,  Mc  +  ml  +  MP,  kg 
=  generalized  coordinates 
=  radius  of  gyration  of  pylon  in  pitch,  roll, 
divided  by  R 

=  radius  of  gyration  of  fuselage  in  pitch,  roll, 
divided  by  R 
=  rotor  radius,  m 

=  beam  mass  divided  by  whole  airframe  mass, 
ml /(ml  +  Me  +  A/p) 

=  transformation  matrix,  Equation  (4) 

=  transformation  matrix,  Equation  (2) 

=  blade  airfoil  velocity  with  respect  to  the  air  mass, 
m/sec 

also  strain  energy,  N-m 
=  velocity  components  of  blade  airfoil  section, 
perpendicular  and  parallel  to  the  chord 
respectively,  m/sec 

=  elastic  deformation  in  x,y,x  directions 
respectively,  m 

=  velocity  vector  of  a  arbitrary  point  on  deformed 
blade,  m/sec 
=  virtual  work,  N-m 

=  distance  along  fuselage,  nose  to  tail,  or 
distance  along  radius  of  rotor,  root  to 
tip,  divided  by  R 


:  rotating  coordinates  fixed  on  the  blade 
=  deformed  coordinates  fixed  on  the  blade 
;  fixed  fuselage  coordinates 
=  dimensionless  fuselage  elastic  degree  of  freedom 
in  vertical  and  lateral  directions 
:  pitch  angle  of  hub,  fuselage,  positive  nose  up, 
rad 

:  roll  angle  of  hub,  fuselage,  positive  advancing 
side  up,  rad 

equilibrium  flapping  angle,  rad 

0o  +  0,  sin^  +  0c  cos  ip 

coning  angle,  rad 

lateral  cyclic  flap  angle,  rad 

longitudinal  cyclic  flap  angle,  rad 

pre-cone  angle 

Lock  number 

first  variation 

scaling  parameter,  (=.l) 

blade  cross-section  principal  axes  coordinates 

vector  of  control  variables 

equilibrium  pitch  angle, 

Bo  +  0,  sin  ip  +  ec  cos  ip  +  Op  (9  -  0pc )  + 
collective  and  cyclic  pitch,  rad 
inflow  ratio 
advance  ratio 

ratio  of  mass  of  pylon  to  mass  of  fuselage, 
MPl(ml  +  Mc) 

ratio  of  lumped  mass  to  the  uniformly 
distributed  mass,  Mp/ml 
air  density,  k g/m3 


=  rotor  solidity 

=  stress  tensor  components,  N/m2 
=  dimensionless  tenssion 

blade  azimuth  angle,  nondimensional  time 
natural  frequency  of  fuselage,  divided  by  0 
frequency  of  ”y”  motion  with  ”x”  boundary 
condition,  divided  by  0;  y  =  z,y,m,L  plunge, 
lateral,  pitch,  roll,  x  =  c,  f  cantilevered, free 
rotor  speed,  rad/sec 
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Figure  1  Fixed,  rotating  and  deformed  Coordinate  systems  of 


blade 


